Introduction à la modélisation statistique bayésienne

Thierry Phénix & Ladislas Nalborczyk
11 juin 2019

Planning

Cours

Cours 1: Introduction à l'inférence bayésienne

Cours 2: Modèle beta-binomial

Cours 3 : Modèle linéaire à 1 prédicteur continu

Cours 4 : Modèle linéaire multivarié

Cours 5: Markov Chain Monte Carlo

Cours 6: Modèle linéaire généralisé

Cours 7: Comparaison de modèles

Cours 8: Prédicteurs catégoriels et Interactions

Cours 9: Modèles multi-niveaux

Cours 10: Data Hackaton

Plan

  • Null Hypothesis Significance Testing
  • Bayes Factor
  • Comparaison de modèles : Overfitting
  • Théorie de l'information
    • Incertitude
    • Divergence
    • Déviance
    • Régularisation
    • Critères d'information
  • Application

Null Hypothesis Significance Testing (NHST)

On s'intéresse aux différences de taille entre hommes et femmes. On va mesurer 100 femmes et 100 hommes.

N = 100
men <- rnorm(N, 175, 10) # 100 tailles d'hommes
women <- rnorm(N, 170, 10) # 100 tailles de femmes
t.test(men, women)

    Welch Two Sample t-test

data:  men and women
t = 2.1475, df = 197.71, p-value = 0.03297
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 0.2655634 6.2364837
sample estimates:
mean of x mean of y 
 173.0157  169.7646 
tibble(
    value = c(men, women),
    gender = c(rep("men", N), rep("women", N))
    ) %>%
    ggplot(aes(x = value, fill = gender) ) +
    geom_histogram() +
    theme_bw(base_size = 20)

plot of chunk unnamed-chunk-3

Null Hypothesis Significance Testing (NHST)

On va simuler des t-valeurs issues de données générées selon l'hypothèse de non-différence entre hommes et femmes.

nSims <- 1e4 # number of simulations
t <- rep(NA, nSims) # initialising an empty vector

t <- replicate(
    nSims, 
    t.test(
        rnorm(N, 170, 10), 
        rnorm(N, 170, 10) 
        )$statistic
    )

La valeur d'un T-test est donnée par la formule :

\[ t = \frac{\mu_W - \mu_M}{S} \]

Où \( S \) est une fonction de la taille des échantillons et des variances des données.

data.frame(t = t) %>%
    ggplot(aes(x = t) ) +
    geom_histogram() +
    theme_bw(base_size = 25)

plot of chunk unnamed-chunk-5

Null Hypothesis Significance Testing (NHST)

La distribution de Student t

Density, distribution function, quantile function and random generation for the t distribution with df degrees of freedom (and optional non-centrality parameter ncp).

Usage

dt(x, df, ncp, log = FALSE)
pt(q, df, ncp, lower.tail = TRUE, log.p = FALSE)
qt(p, df, ncp, lower.tail = TRUE, log.p = FALSE)
rt(n, df, ncp)
df_obs = t.test(men, women)$parameter
data.frame(t = t) %>%
    ggplot(aes(x = t)) +
    geom_histogram(aes(y = stat(density)), 
                   alpha = 0.4) +
    stat_function(fun = dt, 
        args = list(df = df_obs), size = 1.5) +
    theme_bw(base_size = 25)

plot of chunk unnamed-chunk-7

Null Hypothesis Significance Testing (NHST)

alpha <- .05
abs(qt(alpha / 2, df = df_obs) ) # two-sided critical t-value
[1] 1.972035

plot of chunk unnamed-chunk-9

P-values

tobs <- t.test(men, women)$statistic # observed t-value
tobs %>% as.numeric
[1] 2.147452

plot of chunk unnamed-chunk-11

A p-value (in black) gives the probability of observing the data we observed (or more extreme data), given that the null hypothesis is true.

t.test(men, women)$p.value
[1] 0.03297456
tvalue <- abs(t.test(men, women)$statistic)
2 * integrate(dt, tvalue, Inf, df = df_obs)$value
[1] 0.03297456
2 * (1 - pt(abs(t.test(men, women)$statistic),df_obs) )
         t 
0.03297456 

Planning

Cours

Cours 1: Introduction à l'inférence bayésienne

Cours 2: Modèle beta-binomial

Cours 3 : Modèle linéaire à 1 prédicteur continu

Cours 4 : Modèle linéaire multivarié

Cours 5: Markov Chain Monte Carlo

Cours 6: Modèle linéaire généralisé

Cours 7: Comparaison de modèles

Cours 8: Prédicteurs catégoriels et Interactions

Cours 9: Modèles multi-niveaux

Cours 10: Data Hackaton

Plan

  • Null Hypothesis Significance Testing
  • Bayes Factor
  • Comparaison de modèles : Overfitting
  • Théorie de l'information
    • Incertitude
    • Divergence
    • Déviance
    • Régularisation
    • Critères d'information
  • Application

Bayes Factor

Comparaison de deux modèles:

  • \( H_{0}: \mu_{1} = \mu_{2} \rightarrow \delta = 0 \)
  • \( H_{1}: \mu_{1} \neq \mu_{2} \rightarrow \delta \neq 0 \)

\[ \underbrace{\dfrac{p(H_{0}|D)}{p(H_{1}|D)}}_{posterior\ odds} = \underbrace{\dfrac{p(D|H_{0})}{p(D|H_{1})}}_{Bayes\ factor} \times \underbrace{\dfrac{p(H_{0})}{p(H_{1})}}_{prior\ odds} \]

\[ \text{evidence}\ = p(D|H) = \int p(\theta|H) p(D|\theta,H) \text{d}\theta \]

L'évidence en faveur d'un modèle correspond à la probabilité marginale d'un modèle (le dénominateur du théorème de Bayes), c'est à dire à la likelihood moyennée sur le prior… Ce qui fait du Bayes Factor un ratio de likelihood pondéré (par le prior).

Bayes Factor, exemple

On lance une pièce 100 fois et on essaye d'estimer la probabilité \( \theta \) (le biais de la pièce) d'obtenir face. On compare deux modèles qui diffèrent par leur prior sur \( \theta \).

\[ \begin{align} \mathcal{M_{1}} : y_{i} &\sim \mathrm{Binomial}(n, \theta) \\ \theta &\sim \mathrm{Beta}(6, 10) \\ \end{align} \]

\[ \begin{align} \mathcal{M_{2}} : y_{i} &\sim \mathrm{Binomial}(n, \theta) \\ \theta &\sim \mathrm{Beta}(20, 12) \\ \end{align} \]

plot of chunk unnamed-chunk-13

Bayes Factor, exemple


\[ \begin{align} \mathcal{M_{1}} : y_{i} &\sim \mathrm{Binomial}(n, \theta) \\ \theta &\sim \mathrm{Beta}(6, 10) \\ \end{align} \]

\[ \begin{align} \mathcal{M_{2}} : y_{i} &\sim \mathrm{Binomial}(n, \theta) \\ \theta &\sim \mathrm{Beta}(20, 12) \\ \end{align} \]


\[ \text{BF}_{12} = \dfrac{p(D|H_{1})}{p(D|H_{2})} = \dfrac{\int p(\theta|H_{1}) p(D|\theta,H_{1}) \text{d}\theta}{\int p(\theta|H_{2}) p(D|\theta,H_{2}) \text{d}\theta} = \dfrac{\int \mathrm{Binomial}(n, \theta)\mathrm{Beta}(6, 10)\text{d}\theta}{\int \mathrm{Binomial}(n, \theta)\mathrm{Beta}(20, 12) \text{d}\theta} \]

Bayes Factor, exemple

Le Bayes Factor est la nouvelle p-value

Attention à ne pas interpréter le BF comme un posterior odds

Le BF est un facteur qui nous indique de combien nos prior odds doivent changer, au vue des données. Il ne >nous dit pas quelle est l'hypothèse la plus probable, sachant les données (sauf si les prior odds sont 1:1).

Exemple :

  • \( H_{0} \): there is no such thing as precognition
  • \( H_{1} \): precognition exists

On fait une expérience et on calcule un \( BF_{10} = 27 \). Quelle est la probabilité a posteriori de H1 ?

\[ \underbrace{\dfrac{p(H_{1}|D)}{p(H_{0}|D)}}_{posterior\ odds} = \underbrace{\dfrac{27}{1}}_{Bayes\ factor} \times \underbrace{\dfrac{1}{1000}}_{prior\ odds} = \dfrac{27}{1000} = 0.027 \]

Bayes Factor - test d'hypothèse nulle

library(BayesFactor)
ttestBF(men, women)
Bayes factor analysis
--------------
[1] Alt., r=0.707 : 1.312345 ±0%

Against denominator:
  Null, mu1-mu2 = 0 
---
Bayes factor type: BFindepSample, JZS

Prior JZS, en référence à Jeffreys (1961), et Zellner & Siow (1980).

JZS prior

data.frame(x = c(-10, 10) ) %>%
    ggplot(aes(x = x) ) +
    stat_function(
        fun = dcauchy,
        args = list(location = 0, scale = sqrt(2) / 2), size = 1.5) +
    theme_bw(base_size = 20)

plot of chunk unnamed-chunk-15

\[ \mathrm{Cauchy} = \dfrac{1}{\pi \gamma \Big[ 1 + \big(\frac{x-x_{0}}{\gamma}\big)^{2}\Big]} \]

Planning

Cours

Cours 1: Introduction à l'inférence bayésienne

Cours 2: Modèle beta-binomial

Cours 3 : Modèle linéaire à 1 prédicteur continu

Cours 4 : Modèle linéaire multivarié

Cours 5: Markov Chain Monte Carlo

Cours 6: Modèle linéaire généralisé

Cours 7: Comparaison de modèles

Cours 8: Prédicteurs catégoriels et Interactions

Cours 9: Modèles multi-niveaux

Cours 10: Data Hackaton

Plan

  • Null Hypothesis Significance Testing
  • Bayes Factor
  • Comparaison de modèles : Overfitting
  • Théorie de l'information
    • Incertitude
    • Divergence
    • Déviance
    • Régularisation
    • Critères d'information
  • Application

Comparaison de modèles

Deux problèmes récurrents à éviter en modélisation: l'overffitting et l'underfitting.

Comment s'en sortir ?

  • Utiliser des priors pour régulariser, pour contraindre le posterior (i.e., accorder moins de poids à la likelihood)
  • Utiliser des critères d'information (e.g., AIC, WAIC)

Quelle mesure de “fitting” utilise-t-on ?

\[ R^{2} = \frac{\text{var}(\text{outcome}) - \text{var}(\text{residuals})}{\text{var}(\text{outcome})} = 1 - \frac{\text{var}(\text{residuals})}{\text{var}(\text{outcome})} \]

Comparaison de modèles

ppnames <- c("afarensis","africanus","habilis","boisei",
        "rudolfensis","ergaster","sapiens")
brainvolcc <- c(438, 452, 612, 521, 752, 871, 1350)
masskg <- c(37.0, 35.5, 34.5, 41.5, 55.5, 61.0, 53.5)

d <- data.frame(species = ppnames, brain = brainvolcc, mass = masskg)

d %>%
    ggplot(aes(x = mass, y = brain, label = species) ) +
    geom_point() +
    ggrepel::geom_label_repel(hjust = 0, nudge_y = 50, size = 5) +
    theme_bw(base_size = 20) + xlim(30, 70)

plot of chunk unnamed-chunk-16

Comparaison de modèles

Underfitting

\[ \begin{align} v_{i} &\sim \mathrm{Normal}(\mu_{i}, \sigma) \\ \mu_{i} &= \alpha \\ \end{align} \]

mod1.7 <- lm(brain ~ 1, data = d)

plot of chunk unnamed-chunk-18

Comparaison de modèles

Overfitting

mod1.1 <- lm(brain ~ mass, data = d)
(var(d$brain) - var(residuals(mod1.1) ) ) / var(d$brain)
[1] 0.490158
mod1.2 <- lm(brain ~ mass + I(mass^2), data = d)
(var(d$brain) - var(residuals(mod1.2) ) ) / var(d$brain)
[1] 0.5359967
mod1.3 <- lm(brain ~ mass + I(mass^2) + I(mass^3), data = d)
(var(d$brain) - var(residuals(mod1.3) ) ) / var(d$brain)
[1] 0.6797736

Overfitting

plot of chunk unnamed-chunk-21

Comparaison de modèles

Conclusion

  • Augmenter le nombre de paramètres conduit à une amélioration du fitting, \( R^2 \) augmente.

\( \longrightarrow \) Est-ce que l'on cherche vraiment le modèle qui fit le mieux aux données ?

  • Mais augmenter le nombre de paramètres peut conduire à faire des prédictions absurdes !!!

\( \longrightarrow \) Comment choisir le bon nombre de paramètres ?

Il nous faut un outil OBJECTIF de comparaison de modèles !

Planning

Cours

Cours 1: Introduction à l'inférence bayésienne

Cours 2: Modèle beta-binomial

Cours 3 : Modèle linéaire à 1 prédicteur continu

Cours 4 : Modèle linéaire multivarié

Cours 5: Markov Chain Monte Carlo

Cours 6: Modèle linéaire généralisé

Cours 7: Comparaison de modèles

Cours 8: Prédicteurs catégoriels et Interactions

Cours 9: Modèles multi-niveaux

Cours 10: Data Hackaton

Plan

  • Null Hypothesis Significance Testing
  • Bayes Factor
  • Comparaison de modèles : Overfitting
  • Théorie de l'information
    • Incertitude
    • Divergence
    • Déviance
    • Régularisation
    • Critères d'information
  • Applica*tion

Théorie de l'information

  • On voudrait mesurer la distance entre notre modèle et le “vrai” model (i.e., la “réalité”, la nature)…

  • De combien notre incertitude est-elle réduite en apprenant un outcome supplémentaire ?
    Cette réduction est la définition de l'information.

  • Quelle mesure pour l'incertitude ? L'entropie
    S'il existe \( n \) événements possibles, et que chaque événement \( i \) a pour probabilité \( p_{i} \), alors une mesure de l'incertitude est:

\[ H(p) = - \text{E[log}(p_{i})] = - \sum_{i=1}^{n}p_{i} \text{log}(p_{i}) \]

En d'autres termes, l'incertitude contenue dans une distribution de probabilités est la log-probabilité moyenne d'un événement.

Incertitude

Exemple de prédiction météorologique. Si on imagine que la probabilité qu'il pleuve ou qu'il fasse beau (à Lausanne) sont \( p_{1}=0.3 \) et \( p_{2}=0.7 \).

Alors, \( H(p) = - (p_{1} \text{log}(p_{1}) + p_{2} \text{log}(p_{2})) \approx 0.61 \).

p <- c(0.3, 0.7)
- sum(p * log(p) )
[1] 0.6108643

Imaginons que nous habitions à Abu Dabi, et que la probabilité qu'il pleuve ou qu'il fasse beau sont \( p_{1}=0.01 \) et \( p_{2}=0.99 \).

p <- c(0.01, 0.99)
- sum(p * log(p) )
[1] 0.05600153

Entropie croisée

On a donc un moyen de quantifier l'incertitude.
Comment utiliser cette mesure pour quantifier la distance entre notre modèle et la réalité ?

Un modèle avec moins d'incertitude est-il un meilleur modèle ?


Quelle est l'incertitude ajoutée par l'utilisation d'une distribution de probabilités pour décrire… une autre distribution de probabilités ?

Introduisons la notion d'Entropie croisée:

\[ H(p,q) = \sum_{i} p_{i} \log (q_{i}) \]

sum(p * (log(q) ) )

Divergence

  • La Divergence est définie comme l'entropie additionnelle ajoutée en utilisant \( q \) pour décrire \( p \).

\[ \begin{align} D_{KL}(p,q) &= H(p,q) - H(p) \\ &= - \sum_{i} p_{i} \log(q_{i}) - \big( - \sum_{i} p_{i} \log(p_{i}) \big) \\ &= - \sum_{i} p_{i} \big(\log(q_{i}) - \log(p_{i}) \big) \\ \end{align} \]

\[ D_{KL}(p,q) = \sum_{i} p_{i}\big(\text{log}(p_{i}) - \text{log}(q_{i})\big) = \sum_{i} p_{i} \text{log}\bigg(\frac{p_{i}}{q_{i}}\bigg) \]

Divergence

Exemple

Supposons que la “véritable” distribution de nos événements (soleil vs pluie) soit \( p_{1}=0.3 \) et \( p_{2}=0.7 \). Si nous pensons plutôt que ces événements arrivent avec une probabilité \( q_{1}=0.25 \) et \( q_{2}=0.75 \), quelle quantité d'incertitude avons-nous ajoutée ?

p <- c(0.3, 0.7)
q <- c(0.25, 0.75)

sum(p * log(p / q) )
[1] 0.006401457

ATTENTION cette notion n'est pas symétrique

sum(q * log(q / p) )
[1] 0.006164264

Divergence

Courbe de divergence

t <- 
  tibble(
      p_1  = .3,
      p_2  = .7,
      q_1  = seq(from = .01, to = .99, by = .01)) %>% 
  mutate(q_2  = 1 - q_1) %>%
  mutate(d_kl = (p_1 * log(p_1 / q_1)) + (p_2 * log(p_2 / q_2)))
t %>% 
  ggplot(aes(x = q_1, y = d_kl)) +
  geom_vline(xintercept = .3, color = "red", linetype = 2) +
  geom_line(color = "steelblue", size = 1.5) +
  annotate(geom = "text", x = .4, y = 1.5, label = "q = p",
           color = "red", family = "Courier", size = 3.5) +
  labs(x = "q[1]",
       y = "Divergence of q from p")


plot of chunk unnamed-chunk-29

Vers la déviance...

PROBLEME : Nous ne connaissons pas la distribution target (la réalité), à quoi cela peut donc nous servir ?


Astuce: si nous comparons deux modèles, \( q \) et \( r \), pour approximer \( p \), nous allons comparer leurs divergences… Et donc \( \text{E}[\text{log}(p_{i})] \) sera la même quantité pour les deux modèles !

plot of chunk unnamed-chunk-30

Vers la déviance...

On peut donc utiliser \( \text{E}[\text{log}(q_{i})] \) et \( \text{E}[\text{log}(r_{i})] \) comme estimateurs de la distance entre chaque modèle et notre distribution cible. On a donc seulement besoin de la log-probabilité moyenne des modèles. Comme on ne connaît pas la distribution cible, cela veut dire qu'on ne peut pas interpréter cette quantité en termes absolus mais seulement en termes relatifs. Ce qui nous intéresse c'est \( \text{E}[\text{log}(q_{i})] - \text{E}[\text{log}(r_{i})] \).

Rappel : \( D_{KL}(p,q) = H(p,q) - H(p)~~~ \) avec \( H(p,q) = \text{E}[\text{log}(q_{i})] \)

plot of chunk unnamed-chunk-31

Déviance

Pour approximer la valeur de \( \text{E}[ \log(q_{i})] \), on peut utiliser la déviance d'un modèle, qui est une mesure du fit relatif du modèle.

\[ D(q) = -2 \sum_{i} \log(q_{i}) \]

où \( i \) indice chaque observation et \( q_{i} \) est la likelihood de chaque observation.

d$mass.s <- scale(d$mass)
mod1.8 <- lm(brain ~ mass.s, data = d)

-2 * logLik(mod1.8) # compute deviance
'log Lik.' 94.92499 (df=3)

Déviance

# extracting model's coefficients

alpha <- coef(mod1.8)[1]
beta <- coef(mod1.8)[2]

# computing the log-likelihood

ll <- sum(dnorm(
    d$brain,
    mean = alpha + beta * d$mass.s,
    sd = sd(residuals(mod1.8) ),
    log = TRUE)
    )

# computing the deviance

(-2) * ll
[1] 95.00404

In-sample and out-of-sample

La déviance a le même problème que le \( R^{2} \), lorsqu'elle est calculée sur l'échantillon observé. Dans ce cas, on l'appelle déviance in-sample.

Si on est intéressé par les capacités de prédiction de notre modèle, nous pouvons calculer la déviance du modèle sur de nouvelles données… qu'on appellera dans ce cas déviance out-of-sample. Cela revient à se demander si notre modèle est performant pour prédire de nouvelles données.

Imaginons que nous disposions d'un échantillon de taille \( N \), que nous appellerons échantillon d'apprentissage (training). Nous pouvons calculer la déviance du modèle sur cet échantillon (\( D_{train} \) ou \( D_{in} \)). Si nous acquérons ensuite un nouvel échantillon de taille \( N \) issu du même processus de génération de données (que nous appellerons échantillon de test), nous pouvons calculer une déviance sur ce nouvel échantillon, en utilisant les paramètres estimés avec l'échantillon d'entraînement (que nous appellerons \( D_{test} \) ou \( D_{out} \)).

In sample and out of sample deviance

\[ \begin{align} y_{i} &\sim \mathrm{Normal}(\mu_{i},1) \\ \mu_{i} &= (0.15) x_{1,i} - (0.4) x_{2,i} \\ \end{align} \]

On a réalisé ce processus 10.000 fois pour cinq modèles de régression linéaire de complexité croissante. Les points bleus représentent la déviance calculée sur l'échantillon d'apprentissage et les points blancs la déviance calculée sur l'échantillon de test.

Régularisation

Une autre manière de lutter contre l'overfitting est d'utiliser des priors sceptiques qui vont venir ralentir l'apprentissage réalisé sur les données (i.e., accorder plus de poids au prior).

plot of chunk unnamed-chunk-34

Régularisation

Comment décider de la précision du prior ? Est-ce que le prior est “assez” régularisateur ou pas ?

On peut diviser le jeu de données en deux parties (training et test) afin de choisir le prior qui produit la déviance out-of-sample la plus faible. On appelle cette stratégie la cross-validation.

Critères d'information

On mesure ici la différence entre la déviance in-sample (en bleu) et la déviance out-of-sample (en noir). On remarque que la déviance out-of-sample est presque exactement égale à la déviance in-sample, plus deux fois le nombre de paramètres du modèle…

Akaike information criterion

L'AIC fournit une approximation de la déviance out-of-sample:

\[ \text{AIC} = D_{train} + 2p \]

où \( p \) est le nombre de paramètres libres (i.e., à estimer) dans le modèle. L'AIC donne donc une approximation des capacités de prédiction du modèle.

Les limites de l'AIC

l'AIC fonctionne bien uniquement quand :

  • le nombre d'observations \( N \) est largement supérieur au nombre de paramètres \( p \). Dans le cas contraire, on utilise plutôt l'AICc (corrected AIC, voir Burnham & Anderson, 2002; 2004).

  • les priors sont plats (i.e., peu informatifs) ou dépassés par la likelihood (on a beaucoup de données).

  • la distribution a posteriori est une distribution gaussienne multivariée.

Deviance information criterion

Une contrainte de l'AIC est que les priors doivent être plats (i.e., peu informatifs) ou dépassés par la likelihood (on a beaucoup de données). Le DIC est un indice qui ne requiert pas cette condition, en s'accommodant de priors informatifs.

Le DIC est calculé à partir de la distribution a posteriori de la déviance calculée sur l'échantillon d'apprentissage (i.e., \( D_{train} \)).

\[ \text{DIC} = \bar{D} + (\bar{D} - \widehat{D} ) = \widehat{D} + p_{D} \]

où \( \bar{D} \) est la moyenne de la distribution a posteriori \( D \) calculée pour chaque valeur de paramètre échantillonnée, et \( \widehat{D} \) la déviance calculée à la moyenne de la distribution a posteriori. La différence \( \bar{D} - \widehat{D} = p_{D} \) est analogue au nombre de paramètres utilisé dans le calcul de l'AIC.

On ne dira rein de plus ici…

Widely applicable information criterion

Une contrainte d'application de l'AIC et du DIC est que la distribution a posteriori soit une distribution gaussienne multivariée. Le WAIC relache cette condition, en étant souvent plus précis que le DIC.

Un aspect important du WAIC est qu'il est dit pointwise, c'est à dire qu'il considère l'imprécision de prédiction point par point (donnée par donnée), indépendemment pour chaque observation.

On va commencer par calculer la log-pointwise-predictive-density (lppd), définie par:

\[ \text{lppd} = \sum_{i=1}^{N} \log \Pr(y_{i}) \]

où \( \Pr(y_{i}) \) représente la likelihood moyenne de l'observation \( i \) dans le training sample.

En français: la log-densité prédictive point par point est la somme du log de la likelihood moyenne de chaque observation. Il s'agit de l'analogue point par point de la déviance, moyenné sur toute la distribution a posteriori.

Widely applicable information criterion

Application

library(brms)
data(cars)

b <- 
  brm(data = cars, family = gaussian,
      dist ~ 1 + speed,
      prior = c(prior(normal(0, 100), class = Intercept),
                prior(normal(0, 10), class = b),
                prior(uniform(0, 30), class = sigma)),
      iter = 2000, warmup = 1000, chains = 4, cores = 4,
      seed = 6)
 Family: gaussian 
  Links: mu = identity; sigma = identity 
Formula: dist ~ 1 + speed 
   Data: cars (Number of observations: 50) 
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup samples = 4000

Population-Level Effects: 
          Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
Intercept   -17.40      7.07   -30.97    -3.67       1920 1.00
speed         3.92      0.44     3.07     4.75       1783 1.00

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sigma    15.83      1.69    12.98    19.47       2068 1.00

Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample 
is a crude measure of effective sample size, and Rhat is the potential 
scale reduction factor on split chains (at convergence, Rhat = 1).

Widely applicable information criterion

Calcul de \( \text{lppd} = \sum_{i=1}^{N} \log \Pr(y_{i}) \)

Le calcul n'est pas simple

ll <-
  b %>%
  log_lik() %>%
  as_tibble()
dfmean <-
  ll %>%
  exp() %>%
  summarise_all(mean) %>%
  gather(key, means) %>%
  select(means) %>%
  log()

(
  lppd <-
  dfmean %>%
  sum()
)
[1] -206.6836

Widely applicable information criterion

Calcul du nombre de paramètres effectifs

La deuxième partie du calcul du WAIC est le nombre de paramètres effectif, \( p_{WAIC} \). On définit \( V(y_{i}) \) comme la variance de la log-likelihood pour chaque observation \( i \) de l'échantillon d'entraînement.

\[ p_{\text{WAIC}} = \sum_{i=1}^{N} V(y_{i}) \]

dfvar <-
  ll %>%
  summarise_all(var) %>%
  gather(key, vars) %>%
  select(vars) 

pwaic <-
  dfvar %>%
  sum()

pwaic
[1] 3.418859

Calcul de WAIC

Ensuite, le WAIC est défini par:

\[ \text{WAIC} = -2(\text{lppd} - p_{\text{WAIC}}) \]

(WAIC <- -2 * (lppd - pwaic))
[1] 420.205

Le WAIC est aussi un estimateur de la déviance out-of-sample. La fonction waic du package brms permet de faire les calculs.

waic(b)

Computed from 4000 by 50 log-likelihood matrix

          Estimate   SE
elpd_waic   -210.1  6.4
p_waic         3.4  1.2
waic         420.2 12.7

Critères d'information et régularisation

Le DIC et le WAIC peuvent être conceptualisés (au même titre que l'AIC) comme des approximations de la déviance out-of-sample. On remarque que le WAIC produit des approximations plus précises que le DIC, et que l'usage combiné des critères d'information et des priors régularisateurs permet une meilleure approximation de la déviance out-of-sample.

Et après ?

  • Sélection de modèle. On choisit le meilleur modèle un utilisant un des outils présentés et on base nos conclusions sur les paramètres estimés par ce meilleur modèle.

  • Comparaison de modèles. On utilise DIC/WAIC mais aussi les outils de posterior predictive checking discutés précédemment, pour chaque modèle.

  • Moyennage de modèles. On va construire des posterior predictive checks qui exploitent ce qu'on sait des capacités de prédiction de chaque modèle (via le DIC/WAIC).

Planning

Cours

Cours 1: Introduction à l'inférence bayésienne

Cours 2: Modèle beta-binomial

Cours 3 : Modèle linéaire à 1 prédicteur continu

Cours 4 : Modèle linéaire multivarié

Cours 5: Markov Chain Monte Carlo

Cours 6: Modèle linéaire généralisé

Cours 7: Comparaison de modèles

Cours 8: Prédicteurs catégoriels et Interactions

Cours 9: Modèles multi-niveaux

Cours 10: Data Hackaton

Plan

  • Null Hypothesis Significance Testing
  • Bayes Factor
  • Comparaison de modèles : Overfitting
  • Théorie de l'information
    • Incertitude
    • Divergence
    • Déviance
    • Régularisation
    • Critères d'information
  • Application

Comparaison de modèles

Nous allons essayer de prédire les kg par gramme de lait (kcal.per.g) avec les prédicteurs neocortex et le logarithme de mass. Nous allons ensuite fitter 4 modèles qui correspondent aux 4 combinaisons possibles de prédicteurs et les comparer en utilisant le WAIC.

milk <- read.csv("Data/milk.csv")
d <- 
  milk %>%
  drop_na(ends_with("_s"))
rm(milk)

d <-
  d %>%
  mutate(neocortex = neocortex.perc / 100)
head(d)
   X            clade            species kcal.per.g perc.fat perc.protein
1  1    Strepsirrhine     Eulemur fulvus       0.49    16.60        15.42
2  6 New World Monkey Alouatta seniculus       0.47    21.22        23.58
3  7 New World Monkey         A palliata       0.56    29.66        23.46
4  8 New World Monkey       Cebus apella       0.89    53.41        15.80
5 10 New World Monkey         S sciureus       0.92    50.58        22.33
6 11 New World Monkey   Cebuella pygmaea       0.80    41.35        20.85
  perc.lactose mass neocortex.perc neocortex
1        67.98 1.95          55.16    0.5516
2        55.20 5.25          64.54    0.6454
3        46.88 5.37          64.54    0.6454
4        30.79 2.51          67.64    0.6764
5        27.09 0.68          68.85    0.6885
6        37.80 0.12          58.85    0.5885
dim(d)
[1] 17 10

Comparaison de modèles

inits <- list(Intercept = mean(d$kcal.per.g),
              sigma     = sd(d$kcal.per.g))

inits_list <-list(inits, inits, inits, inits)

b6.11 <- 
  brm(data = d, family = gaussian,
      kcal.per.g ~ 1,
      prior = c(prior(uniform(-1000, 1000), class = Intercept),
                prior(uniform(0, 100), class = sigma)),
      iter = 2000, warmup = 1000, chains = 4, cores = 4,
      inits = inits_list,
      seed = 6)

inits <- list(Intercept = mean(d$kcal.per.g),
              neocortex = 0,
              sigma     = sd(d$kcal.per.g))
inits_list <-list(inits, inits, inits, inits)

b6.12 <- 
  brm(data = d, family = gaussian,
      kcal.per.g ~ 1 + neocortex,
      prior = c(prior(uniform(-1000, 1000), class = Intercept),
                prior(uniform(-1000, 1000), class = b),
                prior(uniform(0, 100), class = sigma)),
      iter = 2000, warmup = 1000, chains = 4, cores = 4,
      inits = inits_list,
      seed = 6)
inits <- list(Intercept   = mean(d$kcal.per.g),
              `log(mass)` = 0,
              sigma       = sd(d$kcal.per.g))
inits_list <-list(inits, inits, inits, inits)

b6.13 <-
  update(b6.12, 
         newdata = d,
         formula = kcal.per.g ~ 1 + log(mass),
         inits   = inits_list)

inits <- list(Intercept   = mean(d$kcal.per.g),
              neocortex   = 0,
              `log(mass)` = 0,
              sigma       = sd(d$kcal.per.g))
inits_list <-list(inits, inits, inits, inits)

b6.14 <- 
  update(b6.13, 
         newdata = d,
         formula = kcal.per.g ~ 1 + neocortex + log(mass),
         inits   = inits_list)

Comparaison de modèles

# compute and save the WAIC information for the next three models
b6.11 <- add_criterion(b6.11, "waic")
b6.12 <- add_criterion(b6.12, "waic")
b6.13 <- add_criterion(b6.13, "waic")
b6.14 <- add_criterion(b6.14, "waic")

# compare the WAIC estimates
w <- loo_compare(b6.11, b6.12, b6.13, b6.14,
                 criterion = "waic")

print(w)
      elpd_diff se_diff
b6.14  0.0       0.0   
b6.11 -4.0       2.4   
b6.13 -4.1       1.7   
b6.12 -4.9       2.5   

Pour retrouver les valeurs classiques, il faut faire un petit calcul

cbind(waic_diff = w[, 1] * -2,
      se        = w[, 2] * 2)
      waic_diff       se
b6.14  0.000000 0.000000
b6.11  8.072437 4.894013
b6.13  8.131379 3.465044
b6.12  9.883458 4.952419

Akaike's weights

w[, 7:8] %>% 
  data.frame() %>% 
  rownames_to_column(var = "model_name") %>% 

  ggplot(aes(x    = model_name, 
             y    = waic, 
             ymin = waic - se_waic, 
             ymax = waic + se_waic)) +
  geom_pointrange(shape = 21, color = "steelblue") +
  coord_flip() +
  labs(x = NULL, y = NULL,
       title = "WAIC plot") +
  theme_classic(base_size  = 25) +
  theme(text             = element_text(family = "Courier"),
        axis.ticks.y     = element_blank(),
        panel.background = element_rect())

plot of chunk unnamed-chunk-47

Le poids d'un modèle est une estimation de la probabilité que ce modèle fera les meilleures prédictions possibles sur un nouveau jeu de données, conditionnellement au set de modèles considéré.

\[ w_{i} = \frac{\exp (-\frac{1}{2}\text{dWAIC}_{i})}{\sum_{j=1}^{m} \exp (-\frac{1}{2}\text{dWAIC}_{j})} \]

model_weights(b6.11, b6.12, b6.13, b6.14, 
              weights = "waic") %>% 
  round(digits = 2)
b6.11 b6.12 b6.13 b6.14 
 0.02  0.01  0.02  0.96 

Comparaison des paramètres estimés par nos modèles

library(broom)
my_coef_tab <-
  tibble(model = c("b6.11", "b6.12", "b6.13", "b6.14")) %>% 
  mutate(fit   = purrr::map(model, get)) %>% 
  mutate(tidy  = purrr::map(fit, tidy)) %>% 
  unnest(tidy) %>% 
  filter(term != "lp__")

my_coef_tab %>%
  complete(term = distinct(., term), model) %>%
  select(model, term, estimate) %>%
  mutate(estimate = round(estimate, digits = 2)) %>%
  spread(key = model, value = estimate)
# A tibble: 4 x 5
  term        b6.11 b6.12  b6.13 b6.14
  <chr>       <dbl> <dbl>  <dbl> <dbl>
1 b_Intercept  0.66  0.35   0.71 -1.07
2 b_logmass   NA    NA     -0.03 -0.1 
3 b_neocortex NA     0.46  NA     2.77
4 sigma        0.19  0.19   0.18  0.14

plot of chunk unnamed-chunk-50

Moyennage de modèles

Pourquoi ne conserver uniquement le premier modèle et oublier les autres ? Une autre stratégie consisterait à pondérer les prédictions des modèles par leurs poids respectifs. C'est ce qu'on appelle le moyennage de modèles (model averaging).

  • Calculer le WAIC de chaque modèle
  • Calculer le poids de chaque modèle
  • Simuler des données à partir de chaque modèle
  • Combiner ces valeurs simulées dans un ensemble de prédictions, en utilisant les poids des modèles comme proportions
# we need new data for both the `fitted()` and `pp_average()` functions
nd <- 
  tibble(neocortex = seq(from = .5, to = .8, length.out = 30),
         mass      = 4.5)

# we'll get the `b6.14`-implied trajectory with `fitted()`
f <-
  fitted(b6.14, newdata = nd) %>%
  as_tibble() %>%
  bind_cols(nd)

# the model-average trajectory comes from `pp_average()`
mdl_average <- pp_average(b6.11, b6.12, b6.13, b6.14,
           weights = "waic",
           method  = "fitted",  
           newdata = nd) %>%
    as_tibble() %>%
    bind_cols(nd) 

Moyennage de modèles

Voici les prédictions de tous les modèles considérés, pondérés par leur poids respectif. Comme le modèle b6.14 concentrait quasiment tout le poids, il fait sens que cette prédiction moyennée soit similaire aux prédictions du modèle b6.14.

mdl_average %>%
ggplot(aes(x = neocortex, y = Estimate)) +
  geom_ribbon(aes(ymin = Q2.5, ymax = Q97.5), 
              fill  = "steelblue", alpha = 1/4) +
  geom_line(color   = "darkblue") +
  geom_ribbon(data  = f, aes(ymin = Q2.5, ymax = Q97.5),
              color = "darkgreen", fill = "green", alpha = .1, linetype = 2) +
  geom_line(data = f,
              color = "darkgreen", linetype = 2) +
  geom_point(data = d, aes(y = kcal.per.g), 
             size = 5, color = "orange") +
  labs(y = "kcal.per.g") +
  coord_cartesian(xlim = range(d$neocortex), 
                  ylim = range(d$kcal.per.g))

plot of chunk unnamed-chunk-53

En bleu le modèle le modèle moyenné et en vert le modèle B6.14

Travaux pratiques

Etude des donnés de Howell

  1. Chargement des donnés
    Charger les donnés et créer deux jeux de donnés aléatoirement choisies de même taille.
  2. Crée des modèles polynomiaux de degré 1 à 4
  3. Comparer ces modèles en utilisant le WAIC.
  4. Comparer les rangs des modèles et leurs poids.
  5. Modèle moyen

Pour charger les donnés, utilisez la fonction

read.csv()

Pour selectionner les indice au hasard

sample()

Travaux pratiques

Etude des données de Howell

  1. Chargement des données
  2. Crée des modèles polynomiaux de degré 1 à 4
  3. Comparer ces modèles en utilisant le WAIC.
  4. Comparer les rangs des modèles et leurs poids.
  5. Modèle moyen
Howell1 <- read.csv("Data/Howell1.csv")

d <-
    Howell1 %>%
    mutate(age = scale(age) )

set.seed(1000)
i <- sample(1:nrow(d), size = nrow(d) / 2)

d1 <- d[i, ] # training sample
d2 <- d[-i, ] # test sample

Nous avons maintenant deux dataframes, de 272 lignes chacune. On va utiliser d1 pour fitter nos modèles et d2 pour les évaluer.

Travaux pratiques

Etude des donnés de Howell

  1. Chargement des donnés
  2. Crée des modèles polynomiaux de degré 1 à 4
    Soit \( h_{i} \) les valeurs de taille et \( x_{i} \) les valeurs centrées d'âge, sur la ligne \( i \). Construisez les modèles suivants avec d1, en utilisant brm() et des priors faiblement régularisateurs.
  3. Comparer ces modèles en utilisant le WAIC.
  4. Comparer les rangs des modèles et leurs poids.
  5. Modèle moyen

\[ \begin{align} \mathcal{M_{1}} : h_{i} &\sim \mathrm{Normal}(\mu_{i},\sigma) \\ \mu_{i} &= \alpha + \beta_{1} x_{i} \\ \mathcal{M_{2}} : h_{i} &\sim \mathrm{Normal}(\mu_{i},\sigma) \\ \mu_{i} &= \alpha + \beta_{1} x_{i} + \beta_{2} x_{i}^{2} \\ \mathcal{M_{3}} : h_{i} &\sim \mathrm{Normal}(\mu_{i},\sigma) \\ \mu_{i} &= \alpha + \beta_{1} x_{i} + \beta_{2} x_{i}^{2} + \beta_{3} x_{i}^{3} \\ \mathcal{M_{4}} : h_{i} &\sim \mathrm{Normal}(\mu_{i},\sigma) \\ \mu_{i} &= \alpha + \beta_{1} x_{i} + \beta_{2} x_{i}^{2} + \beta_{3} x_{i}^{3} + \beta_{4} x_{i}^{4} \\ \end{align} \]

Utiliser un intercept initial basé sur la moyenne des hauteurs, et une variance initiale basée sur la variance des hauteurs.

Travaux pratiques

Etude des donnés de Howell

  1. Chargement des donnés
  2. Crée des modèles polynomiaux de degré 1 à 4
    Soit \( h_{i} \) les valeurs de taille et \( x_{i} \) les valeurs centrées d'âge, sur la ligne \( i \). Construisez les modèles suivants avec d1, en utilisant brm() et des priors faiblement régularisateurs.
  3. Comparer ces modèles en utilisant le WAIC.
  4. Comparer les rangs des modèles et leurs poids.
  5. Modèle moyen
inits <- list(Intercept = mean(d$height),
              age       = 0,
              sigma     = sd(d$height))

inits_list <-list(inits, inits, inits, inits)

mod_1 <- brm(
    data = d, family = gaussian,
    height ~ 1 + age,
    prior = c(
        prior(normal(0, 100), class = Intercept),
        prior(normal(0, 10), class = b),
        prior(cauchy(0, 1), class = sigma)),
    iter = 2000, warmup = 1000, chains = 4, cores = 4,
    inits = inits_list)

mod_2 <- brm(
    data = d, family = gaussian,
    height ~ 1 + age + I(age^2),
    prior = c(
        prior(normal(0, 100), class = Intercept),
        prior(normal(0, 10), class = b),
        prior(cauchy(0, 1), class = sigma)),
    iter = 2000, warmup = 1000, chains = 4, cores = 4,
    inits = inits_list
    )

Travaux pratiques

Etude des donnés de Howell

  1. Chargement des donnés
  2. Crée des modèles polynomiaux de degré 1 à 4
    Soit \( h_{i} \) les valeurs de taille et \( x_{i} \) les valeurs centrées d'âge, sur la ligne \( i \). Construisez les modèles suivants avec d1, en utilisant brm() et des priors faiblement régularisateurs.
  3. Comparer ces modèles en utilisant le WAIC.
  4. Comparer les rangs des modèles et leurs poids.
  5. Modèle moyen
mod_3 <- brm(
    data = d, family = gaussian,
    height ~ 1 + age + I(age^2) + I(age^3),
    prior = c(
        prior(normal(0, 100), class = Intercept),
        prior(normal(0, 10), class = b),
        prior(cauchy(0, 1), class = sigma)),
    iter = 2000, warmup = 1000, chains = 4, cores = 4,
    inits = inits_list)

mod_4 <- brm(
    data = d, family = gaussian,
    height ~ 1 + age + I(age^2) + I(age^3) + I(age^4),
    prior = c(
        prior(normal(0, 100), class = Intercept),
        prior(normal(0, 10), class = b),
        prior(cauchy(0, 1), class = sigma)),
    iter = 2000, warmup = 1000, chains = 4, cores = 4,
    inits = inits_list
    )

Travaux pratiques

Etude des donnés de Howell

  1. Chargement des donnés
  2. Crée des modèles polynomiaux de degré 1 à 4
  3. Comparer ces modèles en utilisant le WAIC.
  4. Comparer les rangs des modèles et leurs poids.
  5. Modèle moyen

Utiliser la méthode suivante pour calculer la valeur de WAIC et l'ajouter au modèle

model <- add_criterion(model, "waic")

Utiliser la méthod suivante pour effectuer la comparaison

loo_compare(models, criterion = "waic")

Travaux pratiques

Etude des donnés de Howell

  1. Chargement des donnés
  2. Crée des modèles polynomiaux de degré 1 à 4
  3. Comparer ces modèles en utilisant le WAIC.
  4. Comparer les rangs des modèles et leurs poids.
  5. Modèle moyen
# compute and save the WAIC information for the next three models
mod_1 <- add_criterion(mod_1, "waic")
mod_2 <- add_criterion(mod_2, "waic")
mod_3 <- add_criterion(mod_3, "waic")
mod_4 <- add_criterion(mod_4, "waic")

# compare the WAIC estimates
w <- loo_compare(mod_1, mod_2, mod_3, mod_4,
                 criterion = "waic")

print(w)
      elpd_diff se_diff
mod_4    0.0       0.0 
mod_3  -40.3       9.3 
mod_2 -244.6      19.0 
mod_1 -509.0      21.0 

Travaux pratiques

Etude des donnés de Howell

  1. Chargement des donnés
  2. Crée des modèles polynomiaux de degré 1 à 4
  3. Comparer ces modèles en utilisant le WAIC.
  4. Comparer les rangs des modèles et leurs poids.
  5. Modèle moyen
w[, 7:8] %>% 
  data.frame() %>% 
  rownames_to_column(var = "model_name") %>% 

  ggplot(aes(x = model_name, y = waic, 
             ymin = waic - se_waic, 
             ymax = waic + se_waic)) +
  geom_pointrange(shape = 21, color = "steelblue") +
  coord_flip() +
  labs(x = NULL, y = NULL, title = "WAIC plot") +
  theme_classic(base_size  = 25) +
  theme(text             = element_text(family = "Courier"),
        axis.ticks.y     = element_blank(),
        panel.background = element_rect())

plot of chunk unnamed-chunk-62

Travaux pratiques

Etude des donnés de Howell

  1. Chargement des donnés
  2. Crée des modèles polynomiaux de degré 1 à 4
  3. Comparer ces modèles en utilisant le WAIC.
  4. Comparer les rangs des modèles et leurs poids.
  5. Modèle moyen

Utiliser la méthode suivante pour calculer le poids respectif de chaque modèle

model_weights(models, weights = "waic")

Travaux pratiques

Etude des donnés de Howell

  1. Chargement des donnés
  2. Crée des modèles polynomiaux de degré 1 à 4
  3. Comparer ces modèles en utilisant le WAIC.
  4. Comparer les rangs des modèles et leurs poids.
  5. Modèle moyen

Le poids respectif de chaque modèle

model_weights(mod_1, mod_2, mod_3, mod_4, 
              weights = "waic") %>% 
  round(digits = 2)
mod_1 mod_2 mod_3 mod_4 
    0     0     0     1 

Travaux pratiques

Etude des donnés de Howell

  1. Chargement des donnés
  2. Crée des modèles polynomiaux de degré 1 à 4
  3. Comparer ces modèles en utilisant le WAIC.
  4. Comparer les rangs des modèles et leurs poids.
  5. Modèle moyen
    • Construire les modèles de degré 5 et 6
    • Calculer les poids respectifs de tous les modèles
    • Afficher le modèle moyen
mod_5 <- brm(
    data = d, family = gaussian,
    height ~ 1 + age + I(age^2) + I(age^3) + I(age^4) + I(age^5),
    prior = c(
        prior(normal(0, 100), class = Intercept),
        prior(normal(0, 10), class = b),
        prior(cauchy(0, 1), class = sigma)),
    iter = 2000, warmup = 1000, chains = 4, cores = 4,
    inits = inits_list)

mod_6 <- brm(
    data = d, family = gaussian,
    height ~ 1 + age + I(age^2) + I(age^3) + I(age^4) + I(age^5) + I(age^6),
    prior = c(
        prior(normal(0, 100), class = Intercept),
        prior(normal(0, 10), class = b),
        prior(cauchy(0, 1), class = sigma)),
    iter = 2000, warmup = 1000, chains = 4, cores = 4,
    inits = inits_list
    )

Travaux pratiques

Etude des donnés de Howell

  1. Chargement des donnés
  2. Crée des modèles polynomiaux de degré 1 à 4
  3. Comparer ces modèles en utilisant le WAIC.
  4. Comparer les rangs des modèles et leurs poids.
  5. Modèle moyen
    • Construire les modèles de degré 5 et 6
    • Calculer les poids respectifs de tous les modèles
    • Afficher le modèle moyen
# compute and save the WAIC information for the next three models
mod_5 <- add_criterion(mod_5, "waic")
mod_6 <- add_criterion(mod_6, "waic")

# compare the WAIC estimates
w <- loo_compare(mod_1, mod_2, mod_3, mod_4, mod_5, mod_6,
                 criterion = "waic")
print(w)
      elpd_diff se_diff
mod_6    0.0       0.0 
mod_4   -0.4       1.6 
mod_5   -0.9       1.5 
mod_3  -40.7       9.6 
mod_2 -245.0      18.9 
mod_1 -509.4      20.8 
model_weights(mod_1, mod_2, mod_3, mod_4, mod_5, mod_6, 
              weights = "waic") %>% 
  round(digits = 2)
mod_1 mod_2 mod_3 mod_4 mod_5 mod_6 
 0.00  0.00  0.00  0.32  0.20  0.48 

Travaux pratiques

Etude des donnés de Howell

  1. Chargement des donnés
  2. Crée des modèles polynomiaux de degré 1 à 4
  3. Comparer ces modèles en utilisant le WAIC.
  4. Comparer les rangs des modèles et leurs poids.
  5. Modèle moyen
    • Construire les modèles de degré 5 et 6
    • Calculer les poids respectifs de tous les modèles
    • Afficher le modèle moyen
# we need new data for both the `fitted()` and `pp_average()` functions
nd <- 
  tibble(age = seq(from = -2., to = 3., length.out = 50))

# we'll get the `b6.14`-implied trajectory with `fitted()`
f <-
  fitted(mod_6, newdata = nd) %>%
  as_tibble() %>%
  bind_cols(nd)

# the model-average trajectory comes from `pp_average()`
mdl_average <- pp_average(mod_1, mod_2, mod_3, mod_4, mod_5, mod_6,
           weights = "waic",
           method  = "fitted",  
           newdata = nd) %>%
    as_tibble() %>%
    bind_cols(nd) 

Moyennage de modèles

Voici les prédictions de tous les modèles considérés, pondérés par leur poids respectif. Comme le modèle b6.14 concentrait quasiment tout le poids, il fait sens que cette prédiction moyennée soit similaire aux prédictions du modèle b6.14.

mdl_average %>%
ggplot(aes(x = age, y = Estimate)) +
  geom_ribbon(aes(ymin = Q2.5, ymax = Q97.5), 
              fill  = "steelblue", alpha = 1/4) +
  geom_line(color   = "darkblue") +
  geom_ribbon(data  = f, aes(ymin = Q2.5, ymax = Q97.5),
              color = "darkgreen", fill = "green", alpha = .1, linetype = 2) +
  geom_line(data = f,
              color = "darkgreen", linetype = 2) +
  geom_point(data = d, aes(y = height), 
             size = 2, color = "orange") +
  labs(y = "kcal.per.g") +
  coord_cartesian(xlim = range(d$age), 
                  ylim = range(d$height)) +

plot of chunk unnamed-chunk-69